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ABSTRACT 

We report a new finite-difference scheme for two-dimensional magnetohydrodynamics (MHD) simu- 
lations, with and without rotation, in unstructured grids with quadrilateral cells. The new scheme is 
implemented within the code VULCAN/2D, which already includes radiation-hydrodynamics in vari- 
ous approximations and can be used with arbitrarily moving meshes (ALE). The MHD scheme, which 
' consists of cell-centered magnetic field variables, preserves the nodal finite difference representation 

I of div(B) by construction, and therefore any initially divergence-free field remains divergence-free 

. through the simulation. In this paper, we describe the new scheme in detail and present compar- 

isons of VULCAN/2D results with those of the code ZEUS/2D for several one-dimensional and two- 
[ dimensional test problems. The code now enables two-dimensional simulations of the collapse and 

, explosion of the rotating, magnetic cores of massive stars. Moreover, it can be used to simulate 

["t I ' the very wide variety of astrophysical problems for which multi-D radiation-magnetohydrodynamics 

(RMHD) is relevant. 

Subject headings: magnetohydrodynamics, multi-dimensional radiation hydrodynamics, supernovae 
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1. INTRODUCTION 

\^ ' Magnetic fields play a role, sometimes pivotal, in the evolution and dynamics of astrophysical objects. These include, 
, but are not limited to, molecular clouds, protostars, stellar and planetary dynamos, pulsars, jets from active galactic 

■ nuclei, the solar wind, solar flares, the Earth's and Sun's magnetosheaths, and, importantly, accretion disks of all 

■ kinds. To investigate these environments theoretically has often required multi-dimensional simulation tools that in 
t — , the past have incorporated the relevant physics only to varying degrees. The ideal MHD (magnetohydrodynamic) 

approximation, in which it is assumed the electrical conductivity of the fluid is "inflnite" and, therefore, that the 
magnetic flux is frozen in the matter, is a very good one in many environments in the Universe and obviates the 
necessity to solve separate Boltzmann equations with microscopic couplings for the various charged components. 
■ I Recently, core-collapse supernova explosions and gamma-ray bursts have been added to the list of astrophysical sites 

O ' in which magnetic fields might play an important role. In this paper, we describe a new 2.5-dimensional MHD algorithm 
^ (now incorporated into the radiation hydrodynamic code VULCAN/2D) that maintains the divergence-free character 
^2 of the B-field to machine precision and applies to general structured and unstructured grids. We were motivated in 
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developing this new capability by the intriguing possibility that magnetic fields might play a dynamical role in the 
mechanism of gamma-ray bursts and, possibly, in the explosions of a subset of core-collapse supernovae, as well as by 
an interest in the origin of pulsar and magnetar B-fields. Nevertheless, the new code can be used for a much wider 
set of problems, since along with rotation and magnetic fields in the ideal MHD approximation it also incorporates 
multi-group transport (in both the multi-angle S„ and flux-limited diffusion approximations) and a general gravity 
solver. 

The possibility of hydromagnetic driving of supernova explosions was first explored quantitatively by LeBlanc & 
Wilson (1970). Using numerical simulations that were then the state-of-the-art, they showed that the combination of 
rapid rotation and strong magnetic fields in the pre-coUapse core could lead to the formation of bipolar jets in the post- 
bounce configuration. In that picture, magnetic fields convert rotational energy into the jet energy via the pinch effect. 
The old simulations of LeBlanc & Wilson, however, suffered from crude spatial resolution and did not incorporate 
important developments in the physics of nuclear matter and in the neutrino-matter interaction. They employed 
gray transport and obtained rather weak explosions near '^10^'^ ergs. After this pioneering work, with few exceptions 
(Bisnovatyi-Kogan et al. 1976; Meier et al. 1976; Symbalisty 1984), the possible role of magnetic fields in powering 
or enabling supernova explosions took a back seat to other promising mechanisms. However, recently the subject of 
hydromagnetic interactions during core collapse has been revived by a number of investigators (Akiyama et al. 2003; 
Ardeljan et al. 2005; Yamada & Sawai 2004; Kotake et al. 2004; Ohnishi, Kotake, & Yamada 2005; Akiyama & Wheeler 
2005; Proga 2005; Wilson, Mathews, & Dalhed 2005; Thompson, Quataert, & Burrows 2005; Moiseenko et al. 2006; 
Uzdensky & MacFadyen 2006ab; Masada, Sano, & Shibata 2006), who have investigated the potential of magnetic 
stresses, magnetic buoyancy, magnetic coUimation, magnetic dissipation, and the magnetorotational instability (MRI, 
Balbus & Hawley 1991; Akiyama et al. 2003; Pessah & Psaltis 2005; Pessah, Chan & Psaltis 2006; Obergaulinger, 
Aloy, & Miiller 2006; Etienne, Liu, & Shapiro 2006) to power explosion. In addition, the possibility that magnetic fields 
generated during black hole or neutron star formation might form energetic jets in the context of gamma-ray bursts, 
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X-ray flashes, coUapsars (MacFadyen & Woosley 1999), and hypernovae is receiving increasing attention (e.g., Mizuno 
et al. 2004). Furthermore, the potential role of magnetic winds in powering explosions, or secondary explosions, and 
in spinning down the nascent neutron star or magnetar is now being actively studied (Thompson, Chang, & Quataert 
2004; Bucciantini et al. 2006; Metzger, Thompson, & Quataert 2006). In all these scenarios, rapid rotation and the 
growth of strong magnetic fields due to compression, winding, and/or the MRI during the dynamical phases of core 
collapse play essential roles. Whether progenitors boast the requisite rotation and initial seed fields remains to be 
determined (Heger et al. 2004; Heger, Woosley, & Spruit 2005), but an interesting subset surely must. 

To date, there has been a great deal of MHD code development in astrophysics, mainly focused on simulating 
turbulent magnetic molecular clouds, accretion disks, AGN jets, solar and stellar dynamos, and solar plasmas, but 
lately also on the question of gamma-ray bursts and supernovae. Stone & Norman (1992ab) developed the ZEUS 
suite of codes that employs the method of characteristics (MOC) to properly address MHD waves and incorporates 
the constrained transport (CT) method (Evans & Hawley 1988) to handle the divergence condition on the B-field. 
It includes gray radiation transport using a short-characteristics method. Dai and Woodward (1994) developed an 
MHD Ricmann solver and extended the piecewise parabolic method (PPM) to multidimensional MHD equations. 
They also presented a PPM variant that exactly preserves the conservation laws of magnetohydrodynamics together 
with the divergence-free condition (Dai & Woodward 1998). Gardiner & Stone (2005) have followed ZEUS/MHD 
with a dimensionally-unsplit, 2nd-order, Godunov, grid-based code (Athena) that combines the piecewise parabolic 
method (PPM) for spatial reconstruction and the corner transport upwind (CTU) method of Colella (1990) for time- 
advancing multi-dimensional hydro with the CT method for divergence preservation. Their numerical focus has been on 
making the CTU method compatible with both constrained transport and the finite- volume philosophy of PPM in the 
unsplit context. The result is an accurate and conservative MHD code with great potential to address astrophysical 
problems. Smooth-particle hydrodynamics (SPH) codes have recently been augmented to include magnetic fields, 
using divergence cleaning techniques (Price & Monaghan 2005; Ziegler, Dolag, & Bartelmann 2006) and have been 
employed to study, among other 3D problems, the merger of binary neutron stars (Price & Rosswog 2006). Divergence 
cleaning techniques have been criticized vis a vis staggered-mesh, divergence-free formulations (Balsara & Kim 2004) 
for artificially imposing post-facto cleaning methods to ensure divB — 0, but the approach is still useful for exploratory 
investigations. Aloy et al. (1999ab) have developed a sophisticated special-relativistic, axisymmetric MHD code that 
incorporates as algorithmic elements finite-volume, Godunov, directional splitting, constrained transport, and an HLL- 
type solver. Using a similar code, Obergaulinger et al. (2006) have recently investigated the MRI in the axisymmetric 
core-collapse context. 

There are now several general-relativistic MHD codes that promise to set the standard for astrophysical GRMHD 
simulations in the future. Shibata & Sekiguchi (2005) have constructed a GRMHD code that features fully conservative 
numerical algorithms, dynamic space-times, and the constrained transport method. Noble et al. (2006) have built a 
conservative GRMHD code employing approximate Riemann solvers and a variant of constrained transport. Anninos, 
Fragile, & Salmonson (2005) have produced the GRMHD code, COSMOS-I--I-, that, while it also assumes a fixed 
background space-time, has unstructured grids. COSMOS-I--I- is non-Godunov, uses artificial viscosity, uses finite- 
volume discretization, and has an AMR mode. 

Despite all these recent developments in computational MHD and enormous progress in computing power and in 
numerical techniques, there is still no single code today that can simulate radiation-magnetohydrodynamics (RMHD) 
in full, mainly because a multi-angle, multi-group transport code is expensive, even in 2D. In this paper, we present 
a modest step toward such a code by reporting a new MHD scheme, incorporated into the already operational multi- 
group radiation hydrodynamics code VULCAN/2D (Livne 1993; Livne et al. 2004). Using multi-group, fiux-limited 
diffusion we have obtained and recently published an interesting new explosion mechanism for core-collapse supernovae 
(Burrows et al. 2006ab). We have also explored accretion-induced rotational collapse (Dessart et al. 2006a), convection 
in protoneutron stars (Dessart et al. 2006b), the angle-dependence of the neutrino field and neutrino heating on 
rotation rate (Walder et al. 2005), the mapping between initial and final spins in the core-collapse context (Ott et al. 
2006a), and the gravitational radiation signatures of multi-D supernovae (Ott et al. 2006b). There are several basic 
advantages of VULCAN/2D over other codes (together with some drawbacks) which enable the code to perform in 
difficult situations. One basic feature of VULCAN/2D is its ability to use unstructured grids, which in principle can 
efficiently cover any domain. For example, polar grids become singular at the center of a star and Cartesian grids 
have difficulty resolving the center, while at the same time covering the entire star. In VULCAN/2D, we construct a 
grid by merging an inner curved rectangular with an outer polar grid. The resulting grid resolves the entire star with 
no singularity. In order to maintain this capability we have designed our MHD scheme for unstructured grids, having 
arbitrarily-shaped (but convex) quadrilateral cells. Another approach has been taken by Ardeljan et al. 2000, who 
use a 2D grid composed of triangles only, that also avoids the singularity of polar grids. 

In this paper, we report on the new MHD technique and present a number of test problems that demonstrate 
the ability of the code to simulate 2D magnetohydrodynamics problems with both toroidal and poloidal fields. The 
paper is constructed as follows. In Sj2l we describe in several subsections the basic equations and the finite-difference 
methodology for arbitrary geometry f ij2.ip . the Lagrangean case (i j2.1.ip . and the Eulerian case f ij2.1.2p . We continue 
in fj2T2] with a discussion of the implementation in axially-symmetric geometries fi )2.2p and follow with a discussion 
of the method for evolving the toroidal component with rotation (i i2.2.ip . We end f}2]with a few words on necessary 
modifications to the Courant/CFL condition and boundary conditions ( ^2.2.2^ . Section [3T] presents one-dimensional 
test problems and i )3.2l presents two-dimensional tests, some of them are new. For some of the test problems the 
magnetic stresses are comparable to the pressure stresses. The associated tables and figures support our conclusion 
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(Q that the automaticaUy divergence-conserving MHD scheme we have implemented for general unstructured grids 
in our RMHD code VULCAN/2D provides competitive and reasonably accurate solutions to the dynamic equations 
of MHD in the Newtonian regime. 

2. THE MHD FINITE-DIFFERENCE SCHEME 

In this paper we do not consider gravity, and therefore, omit gravity terms from the equations. We also omit here 
all possible interactions between matter and radiation. This physics, very relevant to the core-collapse problem, is 
already included in VULCAN/2D and has been reported in several papers (Burrows et al. 2006ab; Ott et al. 2006ab, 
Dessart et al. 2006ab) . The MHD equations for incompressible flow in Lagrangean coordinates are (Landau & Lifshitz 
1960; Chandrasekhar 1961): 



dp 
'dt 



pV-v = (1) 

P^ = -Vp+^(VxB)xB (2) 

de d ,1 ^ 

dt=-^dt^-p^ 

|(-) = (-)-Vv. (4) 
dt p p 

Here, e is the specific energy per gram and p{= p{p, e, l^e)) is the pressure, given by the equation of state for a given 
density, p, internal specific energy per gram, e, and electron fraction. Ye (the mean number of electrons per baryon). 
The Lagrangean time derivative ^ can be translated to the laboratory- frame by the transformation ^ — ^ -I- v • V. 

The majority of current numerical schemes for hydrodynamics, including ours, use cell-centered variables for the 
construction of finite-difference approximations to the fluid equations. But unlike Godunov-type schemes, like PPM, 
VULCAN/2D uses staggered differencing where the acceleration and velocity are computed at cell nodes rather than 
at cell centers. In such staggering schemes the acceleration is computed by integrating the pressure times the normal 
over the surface of a control volume around each node. Therefore, we seek a scheme in which the magnetic field is also 
cell-centered, so that the contribution of the magnetic stress to acceleration can be handled in a similar fashion. This 
approach is different, for example, from the choices made in MOC-CT schemes (Stone & Norman 1992b), where the 
magnetic field variables are face-centered (component Bi is face-centered in the i— direction, but cell-centered in the 
perpendicular direction). We do not calculate Lorentz forces directly but rather compute the magnetic stress forces 
from the conservative form of the equations of motion. This is achieved using standard integration schemes, as is done 
for example in Ardeljan (2000). One advantage of our approach is the exact conservation of momentum through the 
simulation. In VULCAN/2D, we use the Cartesian coordinates (x, y) in the planar case and the cylindrical coordinates 
(r, z) in the axisymmetric case. With these choices of coordinate systems, our scheme preserves linear momentum by 
construction, {px,Py in the planar case, Pz in the cylindrical case), and, in the latter case, angular momentum as well. 
We shall discuss the magnetic acceleration in more detail in section 2.2, where the cylindrical case is presented. 

2.1. Differencing the Field Equations: Basic Approach for Arbitrary Geometry 

Our field equation incorporates the flux-freezing principle used in Stone & Norman (1992b). For a given surface, S, 
enclosed by a closed loop, C, we define the fiux: 

j j Bnds. (5) 

If the surface moves with some arbitrary grid velocity the temporal derivative of <I> is given by 

'm 

Here, the time derivative is taken along the trajectory ^ = Vg. We can utilize eq. (6) to evolve a cell-centered 
B-field in both Lagrangean and Eulerian schemes in the following way. Take, for example, the two-dimensional planar 
case, where a cell is some arbitrarily-shaped, convex, quadrilateral in the (x, y) plane. In order to evolve the field 
in time, according to the flux-freezing principle, we need to define two control surfaces in that cell (the number of 
control surfaces must be equal to the number of fleld components) . There are several ways of constructing such control 
surfaces inside a given cell. We considered two such ways - first, connecting opposite mid-edge centers and second, 
connecting diagonals of two opposite nodes. For each such choice we can construct a simple scheme for the evolution 
of the magnetic field. However, the diagonal scheme has been found to be more stable, while at the same time obeying 
the divergence-free condition of the magnetic field in a very elegant way. Therefore, we shall focus on this choice in 
the rest of the paper. 

2.1.1. The Lagrangean Case 

Ignoring for now the toroidal component of the field, B^, the field is described by its two poloidal cell-centered 
components B^^ By. By writing two finite-difference approximations of eq. (6) for the two surfaces one gets two linear 



^ (v - Vg) X Bd/ . (6) 
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equations for the two field components. Let us consider first the Lagrangean case, where Vg = v and the right-hand- 
side of eq. (6) vanishes. Denote the two components of the two surfaces by Ax,i and Aj,,i, z = 1,2. Also, denote 
variables at the beginning of the timestep by superscript n and variables at the end of the timestep by superscript 
n+1. The discretization of eq. (6) in this case yields the following two linear equations for the magnetic field at the 
end of the timestep: 

cn+l Tjn+1 , nn+l r>n+l _ nn r/n , on rtn (ry\ 

s:%'b:+' + = si^b: + si^b; . (8) 

The geometrical factors Sx,i^ Sy^i are the components of the vector nds. 

5y,i = +A^,i. (10) 

It is easy to show that a divergence-free field remains divergence-free if we define a finite diff'erence node-centered 
approximation for the divergence. From the divergence theorem, one gets the approximation 

rfw(B)(P) J J B- nds/V (11) 

where P is a grid node, S is a closed surface surrounding that point and V is the volume contained inside S. For 
our diagonal scheme the surface S is composed of the four diagonals surrounding P (see Fig. 1), and, therefore, the 
integral in eq. (11) is unchanged when the fields in the four relevant cells are evolved by eqs. (7)-(10). 

2.1.2. The Eulerian Case 

In the more general case, where Vg ^ v, the right hand side of eq. (6) represents advection of the B-field across 
the grid and does not vanish. We shall focus on the more common Eulerian case, where Vg = 0. For each of our two 
surfaces we now have to compute the right hand side of eq. (6), and by this obtain two linear equations for the new 
cell-centered components of the field. To do this we need to extrapolate field values from cell centers to the four nodes 
of each cell. Note that for our diagonal scheme the divergence-free condition still holds throughout the evolution, 
because the advection terms cancel out in the integration of eq. (11) over a closed loop. This is true regardless of 
the nodaZ-field values. However, for reasons of stability and accuracy, these nodal values must be carefully defined 
according to monotonicity principles. This is a crucial step where the scheme is different from current MOC-CT 
schemes, which rely on one-dimensional flux-splitting in each direction. Our unstructured grid, which is certainly not 
orthogonal, requires more complex extrapolation in order to ensure accuracy and stability. 

The Eulerian form of our difference scheme for the cell-centered components yields the following two linear equations: 

S:+'b:+^ + S;fB-+^ = SI,B: + SI,b; + M[B^Vy - Byv^l (12) 

S^i'B^' + S-y%'B-+^ = Sl^Bl + Sl^B- + M{B^Vy ByV^], . (13) 

Here, the brackets [BxVy — ByVx]i stand for the difference between the enclosed expression on the two edges of diagonal 
i. The accuracy and stability of the scheme depend upon the exact definition of the nodal values of the magnetic field, 
while the velocity values are already node-centered. We, therefore, have to construct appropriate single-valued nodal 
field components from the cell-centered basic field, and those will be used in the advection terms on the right-hand-side 
of eqs. (12)-(13). The construction consists of the following steps: 

• Compute spatial slopes for the field components according to van Leer's monotonic advection scheme (van Leer 
1979). 

• Extrapolate cell-centered values to the nodes using those slopes. Note that for one- dimensional advection of a 
scalar field this procedure provides a monotonic and second-order-accurate scheme when the advection is done 
upwind. 

• There are now four nodal values for each field component at any internal node p, and each of them is related to 
one diagonal surface in one of the surrounding cells (see diagonals a,b,c,d in Fig. 2). Using the velocity vectors 
at that node define the upwind parameter Oi = cos(v, Ti) for each diagonal, where Ti is a unit vector along 
the diagonal i pointing to the node p. According to the sign of this upwind parameter define the type of each 
surrounding surface to be donor-type if > and acceptor-type otherwise. 

• Usually, there will be two donor-type surfaces (diagonals a and d in Fig. 2) and two acceptor-type surfaces 
(diagonals b and c in Fig. 2) at each (internal) node, but in regions where the grid is non-orthogonal this rule 
may not hold. Without loss of generality assume that surfaces 1, 2 are the two most donor-type surfaces out of 
the four surrounding surfaces, with nodal fields Bi,B2 and normals Ni,N2, respectively. Construct the single- 
valued nodal field vector Bp by solving the two linear equations Bp • Nj = B, • Ni , i=l,2, which imply that the 
normal fluxes on the donor-type surfaces are the same when computed using Bp or the primitive extrapolations 
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The above algorithm deserves a few remarks. First, the choice of donor-type surfaces actuaUy guarantees that the two 
chosen diagonals (surfaces) are not parallel, so that this construction is indeed possible. Secondly, the construction 
of the nodal magnetic vector from the donor-type normal components ensures that the new fluxes are kept bound 
through the advection step. This is essentially the heart of our generalization of van Leer's monotonic upwind scheme 
to the case of magnetic fluxes. Only minor changes are needed for boundary points and we omit discussion here of 
those details. Finally, the same algorithm is applicable for general nonzero grid velocity, which could be used in ALE 
simulations by VULCAN/2D. 

2.2. Implementation in Axially- Symmetric Geometry 

Let us use the common coordinates, cylindrical r,z,(j), for the radial, axial, and azimuthal directions. Assuming 
axial symmetry the azimuthal derivative of any variable vanishes everywhere, but the toroidal components of the 
magnetic field and the azimuthal velocity do not vanish. We first present the equations of motion in a form suitable 
for discretization: 

^ dt "4^^^^ dz dr' 2 dr r ^ dr ^ > 

dv, ^11 d{rB,.B,) 1 d{Bl - - B^^) dp 
^ dt inK dr 2 dz ^ dz ^ ' 

d{rv^) _ 1 I djr^BrB^) d{rB,B^) 
P dt AtiK dr ^ dz ^ ' 

to be compared with Ardeljan et al. (2000) and chapter IX of Chandrasekhar (1961) for the incompressible inviscid 
case. The main reason for writing the magnetic terms in this form is that the equation for the axial momentum and 
the equation for the angular momentum have a conservative form. The difference scheme of VULCAN/2D preserves 
these conservation laws by construction. The impulse and acceleration at each grid point are obtained by integrating 
the above equations over the surface of a local control volume, a cell in the staggered grid, which surrounds the point. 
Thus, momentum conservation laws are preserved exactly by the finite-difference scheme, as the integration provides at 
each surface segment two forces which have equal magnitude and opposite direction. Since this technique is standard in 
hydrodynamical schemes which use staggered grids we omit further details. Note however that the radial momentum 
in cylindrical coordinates is not a conserved quantity. 

Returning to the evolution of the magnetic field we first discuss the poloidal components {B^, B^), while the evolution 
of the toroidal field will be described separately in §2.2.11 The numerical approximation of those equations is obtained 
from the flux-freezing principle in a way similar to the planar case. We first take a cell in the (r, z) plane with the two 
surfaces defined by the diagonals connecting opposite nodes. The cell represents a torus around the axis of symmetry 
and the surfaces are actually closed ring-like sheets centered around the axis of symmetry. In the Lagrangean case, 
the resulting two linear equations for {Br, B^) are then 

on+l on+l I c"i+l D"+l C" on i C" (^ Q \ 

where the geometrical factors are 



Sr,, = -27rr,A^,, (19) 

Sz,^^+2^:r,l\r,^, (20) 

ri being the radius at the center of each surface. We see that the only difference with the planar case is the factor 
of 27rr multiplying each geometrical factor. Likewise, in the Eulerian case we get equations similar to eqs. (12)-(13) 
with appropriate 27rr factors. Therefore, the numerical technique for evolving the poloidal field is almost identical to 
that in the planar case. 

2.2.1. The Toroidal Component and Rotation 

For the time evolution of the toroidal component we need an extra equation. In the pure Eulerian case, which is our 
most common case, the equation for the toroidal field takes the form 

dH d d 

= jriBrV^ - B^^r) + T-{BzV^ - B^Vz) . (21) 

dt dr dz 

Using flux freezing again, we define an extra surface which extracts the toroidal field. This is just the rectangular image 
of the cell in the (r, z) plane. Integrating eq. (21) over this surface yields the desired equation. After integration, the 
right-hand-side gives edge-centered advection terms which we compute using the standard upwind van Leer scheme. 



2.2.2. Boundary Conditions and Time Steps for MHD Problems 
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We employ an explicit version of VULCAN/2D hydro and, thus, must restrict the timestep according to stability 
constraints. We use the CFL condition (Courant & Friedrichs 1948) updated to account for the propagation of Alfvcn 
and fast magnetosonic waves. In practice, we use 

At < min (^^) , (22) 

where |ti| is the local velocity (accounting for any rotational component as well) and C/ is the fast-magnetosonic 
speed (C| — v\ + C'^), where va is the Alfven speed and Cs the sound speed. Note that in core-collapse supernova 
simulations, incorporating modest initial magnetic fields does not reduce noticeably the time increment, which after 
bounce is of the order of a microsecond. VULCAN/2D can also operate in an implicit mode where the Courant 
condition for acoustic waves is removed. The MHD solver however is explicit in nature. This implies that, depending 
on the field strength, the timestep in future implicit simulations could be restricted by the speed of Alfven waves. 

As to boundary conditions, they can be rather complicated if one wishes to simulate general MHD problems in 
arbitrary geometries. Since we are interested mainly in the core-collapse problem, and problems of a similar nature 
around compact objects, we choose to sacrifice the general case in favor of the most simple treatment. In our current 
implementation, the dynamical region in embedded in a large enough ambient region where either practically nothing 
changes during the simulations, or the outer boundary is not allowed to move. Therefore, under those restricted 
assumptions, a divergence- free initial field remain stationary in the outer regions of the flow. 

3. TEST PROBLEMS 

A number of one-dimensional and two-dimensional test problems have been chosen to verify the code, in both planar 
and cylindrical geometries. We have computed some of the test problems suggested in Stone et al. (1992) and in 
Stone & Norman (1992b). We have also designed a few new test problems which are relevant to dynamical collapse. 
The results of each VULCAN/2D simulation are compared to a simulation done under identical conditions using the 
pubhcly available MHD code ZEUS/2D (Stone et al. 1992; Stone & Norman 1992b). 

3.1. One-Dimensional Test Problems 

In this section, we present one-dimensional planar and cylindrical test problems. For this purpose, we constructed a 
one-dimensional Lagrangean magnetohydrodynamic version, VULCAN/ID, based upon the MHD formalism outlined 
in this paper and upon the earlier work of Bisnovatyi-Kogan et al. (1976). VULCAN/ID is a variant that assumes that 
all the hydrodynamical and magnetic variables are functions of (r, t) alone. We use this Lagrangean code for testing 
our basic flux freezing approach in constructing the finite-difference scheme, and also for post comparisons between 
VULCAN/2D and VULCAN/ID. In most cases, the results of our VULCAN/2D (Eulerian) simulations are compared 
either to similar results obtained using VULCAN/ID (Lagrangean) or to those obtained using ZEUS/2D (Eulerian). 

3.1.1. Advection of a Magnetic Pulse 

This problem addresses the one-dimensional advection of a pulse of transverse magnetic field in Cartesian geometry, 
following the procedure presented in Evans & Hawley (1988) and in Stone et al. (1992). The velocity is fixed and 
uniform throughout the grid and the exact solution is just a translation of the initial signal. The test probes the 
magnitude of the numerical dispersion and diffusion of the advection scheme employed in the MHD solver, that should 
operate exactly like the same advection scheme of a scalar field. 

Our initial conditions are a one-dimensional grid of 500 zones, equally spaced to tile the x-direction from to 100. 
We use a Courant number of 0.5 and the equation of state is that of an ideal gas. Density and pressure are set to a fixed 
value everywhere, a value which does not matter in practice since we perform pure advection.. The velocity is unity 
everywhere, and we start with a non-zero transverse magnetic field between a; = 10 and x = 20. Boundary conditions 
are periodic. We run the simulation until this pulse has advected five times its initial width. We show in Fig. [3] the 
results for the ZEUS/2D simulations using the donor-cell scheme and the van Leer scheme and for simulations using 
VULCAN/2D. For comparison, we overplot the initial (square) pulse shifted along the x-direction by vt. VULCAN/2D 
results are very close to ZEUS/2D with the van Leer scheme (the two curves, essentially, overlapping). Given this 
good agreement, we do not investigate the distribution of the current density nor estimate errors as shown in Figs. 2-3 
in Stone et al. (1992). As expected, van Leer's scheme introduces a dramatic improvement of the results compared to 
donor cell. Another higher-order variant, the PPA scheme of ZEUS/2D, improves only a little on van Leer's scheme, 
while introducing non-monotonic unphysical oscillations (not shown here). 

3.1.2. Brio- Wu Test Problem 

The Brio-Wu test problem (Brio & Wu 1988) is a one-dimensional, planar, MHD Ricmann problem. Analytic 
solutions to MHD Riemann problems are discussed in Ryu & Jones (1995). The particular Brio-Wu problem belongs 
to a subclass of compound states, where the analytic solution is ambiguous. We, therefore, compare our results to the 
numerical solution of ZEUS/2D. The 2:-range extends from to 100, and the initial parameters on the two sides of the 
initial discontinuity at x = 50 are: 



fx < 50 : p = 1.000, P = 1.0, = 0.75, By = +1 
\a; > 50 : p = 0.125, P = 0.1, B^ = 0.75, By = -1 
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The equation of state is of an ideal gas with ^ — 2. Boundary conditions are reflecting. To illustrate the behavior 
of VULCAN/2D vis d vis ZEUS/2D, we follow the same procedure found in Stone et al. (1992) and perform tests 
at different resolutions, assessing qualitatively and quantitatively the errors and how these vary with resolution. The 
Courant number for all simulations is set to 0.5. 

First, in Fig. [H we present the solution of the Brio-Wu problem at t = 10 computed with VULCAN/2D (black) 
and ZEUS/2D with the van Leer scheme (red). The agreement seen in the graphs is excellent. For a quantitative 
comparison of ZEUS/2D with the van Leer scheme, VULCAN/2D, and the analytic Brio-Wu results, we log in Tabled] 
the values predicted by each for the density, p, the pressure, P, the velocities, 14 and Vy, and the B-field component, 
By. We provide these values at strategic locations in the fluid, which, from left to right, are the fast rarefaction (FR; 
X = 45), the shock discontinuity (SC; x = 48), the left and right sides of the contact discontinuity (CDi and CDj-; 
X — 55), the fast rarefaction (FR; x = 80), as well as the left and right initial boundary values for reference. 

Following Stone et al. (1992), using self-convergence testing, we now estimate errors and convergence rates for both 
ZEUS/2D and VULCAN/2D in this specific Brio-Wu test. The simulations presented in Fig.|4]were done with TV = 512 
zones, uniformly spaced in the x-direction. We performed a series of similar calculations with different resolutions: 
N =64, 128, 256, 1024, and 2048. We then calculated an Li error norm defined by 



Li = 



J I /027V (a;) - Pn{x) | dx 



Jdx 



(23) 



where P2n{x) {pn{x)) represents the density in the run with 2N (N) zones, both taken at t = 10 and at the same 
location x. Note that here Li is a density error norm and that it has the units of a density, rather than being 
dimensionless. Using this measure, we seek the rate at which the error varies with resolution, rather than the absolute 
level of that error. Having six levels of resolution, we compute Li for five consecutive pairs and show the results in 
the left panel of Fig. [5l The convergence rate is similar for both simulations, i.e., ZEUS/2D with van Leer's scheme 
and VULCAN/2D. 

Finally, we display in the right panel of Fig.[5]the relative density errors for the simulations at iV = 256 (broken) and 
TV = 1024 (soUd) zones for VULCAN/2D (Black) and ZEUS/2D/van Leer (red). Stone et al. (1992) noted that such 
an MHD shock-tube test, with a sharp discontinuity, forces code accuracy to be first-order. Our comparisons between 
codes that are completely different in design provide an instructive look into code performance and the accuracy of 
the results generated. 

3.1.3. Propagation of an Alfven Wave in ID 

Here, we examine the propagation of a ID Alfven wave, defined as a transverse perturbation on a given constant 
field Bx = Bq, embedded in a homogeneous medium of constant density p. The solution of such a wave as a function 
of x is naturally (Landau & Lifshitz 1960): 

Vy{x, t) = VQCos{kx — ujt) (24) 
By{x,t) ^ -^/Iirpvy . (25) 



The perturbations Vy and By propagate along the x direction with the Alfven speed va = ■ We evolve the 

solution at t — with fc = 2 in the interval [0, 27r] for five periods and compare the profiles at each integer period to the 
initial configuration. Table [2] lists the LI error of By for various resolutions (TV) and 5 integer periods. As functions 
of resolution, our errors behave similarly to those presented in Table H of Toth (2000) for the traveling wave, and 
show a second-order rate of convergence. As a function of time, we observe a linearly growing error in each column 
in Table [5] due to the fact that our scheme, like most other schemes, is only first-order accurate in time. Moreover, 
in more realistic problems which involve advection and/or acoustic waves, second-order accuracy usually degrades to 
nearly first-order. This can be seen in Toth's Table II for the standing wave and in his Tables V-VIII, where the 7 
tested schemes show a roughly first-order rate of convergence with spatial resolution. 

3.1.4. Transport of Angular Momentum in a Rotating Disk 

This one-dimensional problem in cylindrical geometry tests the evolution of a toroidal field, together with the 
transport of angular momentum in a differentially rotating disk. Here, we start with differential rotation in the radial 
direction and a corresponding radial magnetic field. This configuration is more typical of astrophysical problems than 
is employed in the Braking of aligned rotators problem of Stone et al. (1992). The stresses of the radial magnetic field 
transport angular momentum from the inner parts outward, evolving towards uniform rotation, while a toroidal field 
grows from zero to significant magnitude. Bisnovatyi-Kogan et al.(1976) suggested this mechanism for an MHD-driven 
supernova, but they need extremely strong initial fields to generate an explosion. 

We first set an initial disk configuration in rotational equilibrium, where the centrifugal acceleration is balanced by 
pressure only, namely, 

-^^rcj'ir). (26) 
p or 
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Using an isentropic ideal gas with adiabatic index 7, p = Ap^ , and a rotation law of the form 

-W = ^, (27) 

one can integrate the hydrostatic equation to get the entire disk configuration as a function of the radius. We computed 
this test problem, with the parameters p{r = 0,t = 0) = 100, /9(r = 0, < = 0) = 1,7 = 5/3, ojq = 10 (in dimensionless 
units), using our one dimensional Lagrangean code VULCAN/ID and VULCAN/2D in the Eulerian mode. The radius 
of the disk is Vmax = 10 and the width of the initial rotation profile is L = rmax/4:. In both simulations, we have used 
500 radial zones and fixed, rigid, outer boundary conditions. Without a magnetic field the disk remains static and 
stable in both simulations, showing that the schemes do not introduce spurious numerical errors. For the real test we 
turn on a radial magnetic field of magnitude rBr = 100 and evolve the configuration until t=0.5. Figure [5] displays 
the toroidal field (left) and the angular velocity at three times (0.0,0.3,0.5). The pluses, showing the results of 
VULCAN/2D, are plotted for only one fifth (100) of the zones in order not to hide the solid lines. The excellent 
agreement between the two codes lends credence to the accuracy and stability of our 2D Eulerian scheme, which is 
much more complicated than its ID Lagrangean version. 

3.2. Two- Dimensional Test Problems 

There are only a few published dynamical test problems which are two-dimensional in nature. Therefore, we 
constructed a number of such problems, which are close to astrophysical problems in their dynamics, and compared 
the results of VULCAN/2D with those of ZEUS/2D. We focus on problems which involve supersonic flow and strong 
shock waves, rather than on steady-state problems like wind solutions. In the absence of analytic solutions code-to- 
code comparison is the only way to validate multi-dimensional simulations. However, such comparisons should be 
done with great care. In a few previous attempts, it was found that significant differences between codes persist even 
when high-order schemes are used with maximum available resolution. Examples are the "Santa Barbara Cluster 
Comparison Project" (Frenk et al. 1999), the "Richtmyer-Meshkov Instability Study" of Holmes et al. (1999) and the 
"Rayleigh- Taylor Instability Study" of the "Alpha-Group" (Dimonte et al. 2004). The differences between codes arise 
from the tendency of unstable flows to form smaller and smaller structures due to vorticity generation. Therefore, 
one can anticipate finding good agreement between two codes only on resolved structures, for which there are at least 
10-12 grid zones. In other words, there should be no point-wise convergence in simulations of unstable flows in regions 
with small-scale fragmentation and shear (see also discussions by Calder et al. 2002). Consequently, self-convergence 
tests in the sense described in §3.1.21 are not a common practice in multi-dimensional simulation. 

3.2.1. Exploding Sphere/Blast Wave 

This axisymmetric test problem consists of the expansion of a high pressure sphere embedded inside an ambient 
cloud with initially uniform magnetic field. The specific setup used is (in normalized dimensionless units): 



r 0.01 < i? < 0.2 : p= 100, P = 10^ Br ^0,B^ = 140 
|o.20 < i? < 1.0 : p= 1.0, P ^ 0.60, Br = 0, = 140 , 

where R is the spherical radius, and Br and B^ refer to the cylindrical components of the magnetic field (spherical 
geometry is, however, used for the simulation). The problem was evolved with an ideal gas equation of state with 
7 = 5/3 and the Courant number used is again 0.5. The spherical grid is logarithmically spaced in radius with nr 
zones between 0.01 and 1.0 and uniformly spaced in angle with nt zones between and tt. Boundary conditions are 
reflecting and the simulation is stopped before the blast wave hits the outer boundary. 

We show in Fig. [7] a montage of the density (top- left panel), pressure (top-right panel), and the cylindrical r and z 
components of the magnetic field (bottom-left and bottom-right panels, respectively) at the final time t = 0.0035 in a 
high-resolution simulation with nr — 280 and nt ~ 240. Each panel is broken into two halves, with the VULCAN/2D 
results on the right and the ZEUS/2D results on the left. To give a better rendering of the dynamics and of the 
poloidal field morphology, we overplot velocity vectors (unsaturated; maximum length corresponds to a value of 165 
or a physical length of 10% of the width of the display). Magnetic field lines are drawn starting from z — —0.95 r^^^x, 
equally spaced every 0.1r,„ax in the horizontal direction. The different paths followed by these lines between the two 
sides are indicative of differences in the magnetic field computed by VULCAN/2D and ZEUS/2D. Overall, lines bunch 
up in the shock region where they accumulate due to flux freezing with the mass. 

At the qualitative level, the agreement between the two codes is excellent. The colormaps of the components of the 
poloidal field show the compression of the vertical field Bz and the generation of the radial component due to fiux 
freezing. Because the flow fleld of this problem is irrotational, we performed a resolution study in the sense described 
at the previous section. Since in this test problem there are only large-scale structures without shear we expect to see 
some kind of convergence. We show in Table [3] the results at the final time for the extrema of the density, pressure, 
Vr, Vz, Br, and B^- Figure [7| shows results for the highest-resolution simulation, while in Table [3] we give additional 
data for low resolution (nr = 71 and nt = 61) and medium-resolution [nr = 141 and nt = 121) tests. All quantities 
converge at higher resolution where they differ between the two codes by no more than a few percent. However, 
at low resolution, differences can be of the order of 20-30%. Note, however, that in general the low-resolution and 
medium-resolution values obtained by VULCAN/2D are closer to their high-resolution values than those obtained by 
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ZEUS/2D. Using the highest resolution model as the reference model, the Li error norm for the density (see discussion 
in the Brio-Wu test above) gives a value for the low- and medium-resolution simulations, respectively, of 0.42 and 0.22 
for ZEUS/2D, 0.24 and 0.19 for VULCAN/2D. OveraU, ZEUS/2D and VULCAN/2D compare well for this explosion 
test, both qualitatively and quantitatively. 

3.2.2. Imploding Sphere 

In this test, the initial spherically-symmetric configuration is split between an outer low-density, high-pressure shell 
and an inner high-density, low-pressure region. The entire domain is threaded by a uniform vertical field Bz- The 
initial setup is given by (in normalized dimensionless units): 

r 0.04 < i? < 0.8 : p = 10, P = 1.00, Br = 0, = 40 
\0.80 < i? < 2.0 : p= 1.0, P = 1000, Br = 0, = 40 . 

The problem was evolved with an ideal gas equation of state with 7 = 5/3, and in the highest resolution version, 
we employ nr = 360 logarithmically-spaced radial zones and nt — 240 uniformly-spaced angular zones on a full tt 
sector. The initial configuration evolves so that the outer high-pressure region implodes, compressing the inner region 
right down to the inner boundary at R = 0.04, which then reexpands. The interaction of the magnetic field with 
the implosion leads to fragmentation of the shock wave and to the formation of many small-scale structures. The 
ratio of the gas pressure to the magnetic pressure (/3) reaches unity inside near bounce. Note that for such tests with 
ZEUS/2D a large latitudinal flow developed at mid-latitude in the inner-boundary "shell" promptly after the start of 
the simulation. We cured this anomaly in this otherwise static inner region by modifying the base boundary condition 
in ZEUS/2D to enforce a null latitudinal velocity in the ghost zones and in the first active zone. 

In a first test, we computed the implosion without rotation, for which we show in Fig. [8] results for the density (top 
left), pressure (top right), Br (bottom left) and B^ (bottom right) at time t = 0.05 and in Fig. [9] at time t = 0.07. 
The shock bounces at the inner rigid surface shortly after t — 0.05, and we evolve the flow until t = 0.1. Even before 
shock bounce, a complex pattern is formed due to the interaction between the vertical magnetic field and the radially 
converging flow (Fig. [5]). After bounce, the flow fragments on smaller scales and forms complex structures by non-radial 
shock interactions (Fig. [9]). The two codes agree well on scales which are resolved by more than a few zones. Although 
VULCAN/2D is not designed to conserve energy exactly (with or without magnetic fields) the total energy in this 
simulation is conserved to an accuracy of ~0.2 percent; the magnetic energy accounts for roughly 5 percent of the 
total. Both codes preserve accurately the left-right symmetry of Bz and the anti-symmetry of Br- Those symmetries 
are seen in Fig. [8] and subsequent figures, but can best be rendered by the morphology of the poloidal field lines (top 
right panels of Figs. [8] and [9] and in the right panels of Figs. [TOl and fTTj) . 

VULCAN/2D seems more dissipative than ZEUS/2D on the smallest scales and at sharp discontinuities. The slight 
differences at these smallest scales and at interfaces could be consequences of the fact that VULCAN/2D does not 
currently use any "sharpening" techniques on discontinuities. The agreement between the codes persists all the way 
to the end of the simulation at t = 0.1 (not shown here). In Table 21 we present the integrated energy components 
at several times. We can see that the total energy is conserved at the same level of accuracy in the two codes. The 
agreement between separate components of the energy budget is also good, where the largest difference appears in the 
kinetic energy. Since this difference appears already at early times, and then remains roughly constant, we speculate 
that the origin of this difference is a different treatment of the initial shock wave that emerges from the large initial 
pressure discontinuity. 

A second similar test uses the same initial input, but with initial rotation. The initial rotation speed is given by 
v^(r) = -^j^rjA)'^ choose = 10 and A = 0.5 (r is the cylindrical radius). In Fig. [TOl we present comparisons 

of the angular velocity O (left column) and the toroidal magnetic field B^ (right column) before bounce [t = 0.05; 
top row) and after bounce {t = 0.07; bottom row). The initial field in this simulation is purely poloidal, so that the 
non-zero toroidal field stems from the winding of the poloidal component through differential rotation. The complex 
2D configuration of the angular velocity leads to both positive and negative toroidal components of the magnetic field, 
depending on location (height above the inner boundary and hemisphere). Again, the agreement between the codes is 
good for all structures resolved by more than a few zones. In Table [51 we present the integrated energy components at 
the times used for the non-rotating case. The conservation of total energy and the agreement between different energy 
components are of the same quality as in the previous case. Note, however, that while VULCAN/2D preserves the 
total angular momentum exactly ZEUS/2D does not, though the deviation is only a few percent. 

All simulations presented so far were set up with intermediate strength initial B-fields, and, indeed, comparable 
magnetic and gas pressures were obtained only transiently in the simulations, e.g., during the bounce epoch in the 
implosion test. To test the high B-field case, we repeated the high-resolution non-rotating implosion test discussed 
earlier, but this time increased the initial magnetic field magnitude from 40 to 120 and 300. The magnetic field 
pressure then dominates the gas pressure. We present in Fig.[TT]a montage of colormaps of the density (left) and the 
quantity /3piasma ~ P/{B'^/8t^) (right), the ratio of thermal and magnetic pressures, for the cases with initial fields 
Bz = 120 (top; t = 0.06) and Bz — 300 (bottom; t = 0.04). Note how the material is more and more confined to 
motion along field lines as the field strength increases, preventing any perpendicular motion in the highest i?^ -value 
case. Here, as before, differences are noticeable between VULCAN/2D and ZEUS/2D, for example in the extrema 
reached by the density and /3piasma (or the gas pressure and the magnetic field components) , but overall the agreement 
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is quite satisfactory. Simulations at high magnetic field strength imply large Alfven speed and, thus, reduced timestep. 
Because such simulations are costly, we stopped those simulations earlier than the final times of the intermediate 
strength cases. 

4. CONCLUSIONS 

We have developed an ideal MHD scheme for general unstructured grids in two dimensions, plus rotation. The 
scheme consists of flux freezing of cell-centered magnetic vectors in both Lagrangean and Eulerian realizations. For 
the Eulerian mode, we generalized the van Leer second-order advection scheme for cell-centered fluxes in a unique way. 
The scheme is now incorporated into the (Newtonian) radiation hydrodynamics code VULCAN/2D and can be used 
for future two-dimensional RMHD simulations of rotating core-collapse supernovae and the early phases of gamma-ray 
bursts, as well as for simulations of a variety of other multi-D RMHD astrophysical problems. 

Via a number of one-dimensional and two-dimensional test problems, we find that the accuracy of the scheme is 
comparable to that of high-order MOC-CT schemes used, for example, in the code ZEUS/2D. It is, however, impossible 
to prove a formal order of accuracy for our scheme because it does not consist of one-dimensional sweeps. In actual 
test problems we see that VULCAN/2D is somewhat more dissipative on small, unresolved, scales than ZEUS/2D, 
but the agreement between the two codes on structures that are resolved by more than a few zones is good. The 
extra dissipation of VULCAN/2D is probably related to numerical details that are not related to the MHD scheme, 
and do not affect resolved scales. The MHD scheme is designed to conserve the nodal finite-difference representation 
of div(B) automatically to machine accuracy at all times. Moreover, it conserves linear and angular momenta in the 
axisymmetric cylindrical case. 
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TABLE 1 

Values of selected quantities for the Brio-Wu test at t = 10. 



Variable 


Source 


Left 


FR 


SC 


CDi 




FR 


Right 


P 


VULCAN/2D 
ZEUS/2D 
BW 


1.0 
1.0 


0.664 
0.660 
0.676 


0.839 
0.840 


0.704 
0.699 
0.697 


0.242 
0.250 


0.116 
0.116 


0.125 
0.125 


P 


VULCAN/2D 
ZEUS/2D 
BW 


1.0 
1.0 


0.441 
0.435 
0.457 


0.719 
0.712 


0.505 
0.501 
0.516 


0.505 
0.501 


0.086 
0.086 


0.1 

0.1 




VULCAN/2D 
ZEUS/2D 
BW 


0.0 
0.0 


0.663 
0.672 
0.637 


0.485 
0.488 


0.596 
0.603 
0.599 


0.596 
0.603 


-0.264 
-0.279 


0.0 
0.0 


Vy 


VULCAN/2D 
ZEUS/2D 
BW 


0.0 
0.0 


-0.247 
-0.251 
-0.233 


-1.19 
-1.16 


-1.58 
-1.58 
-1.58 


-1.58 
-1.58 


-0.185 
-0.196 


0.0 
0.0 


By 


VULCAN/2D 
ZEUS/2D 
BW 


1.0 
1.0 


0.568 
0.562 
0.585 


-0.190 
-0.195 


-0.537 
-0.539 
-0.534 


-0.537 
-0.539 


-0.892 
-0.886 


-1.0 
-1.0 


Note. — Hydromagnetic variables at selected locations in the BW solution, which, from left to right, are the fast 
rarefaction (FR; x — 45), the shock discontinuity (SC; x — 48), the left and right sides of the contact discontinuity (CDi 
and CDr; x — 55), the fast rarefaction (FR; x — SO), as well as the left and right initial boundary values for reference (see 
Fig^. Note how these values closely match. 




TABLE 2 

LI ERRORS IN By FOR THE AlFVEN WAVE TEST 






Periods 


N=16 






N=32 




N=64 




N=128 


P=l 


0.0086 






0.0021 




0.00054 




0.00015 


P=2 


0.0171 






0.0041 




0.00106 




0.00028 


P=3 


0.0256 






0.0062 




0.00158 




0.00041 


P=4 


0.0339 






0.0083 




0.00210 




0.00054 


P=5 


0.0423 






0.0103 




0.00262 




0.00068 



Note. — The value of N indicates the number of zones per wavelength 



TABLE 3 

Values of selected quantities for the explosion test at t = 0.0035, using low, medium, and high 

RESOLFTIOX. 



V <Xi lOiDiC 


Code 


Low resolution 


Medium resolution 


High resolution 






Mill. 


Max. 


Min. 


Max. 


Min. 


Max. 


p 


VULCAN/2D 


1.00 


3.53 


1.0 


4.48 


1.0 


5.4 




ZEUS/2D 


1.02 


2.83 


0.99 


4.12 


0.97 


5.3 


P 


VULCAN/2D 


0.60 


10157 


0.6 


10119 


0.6 


10525 




ZEUS/2D 


2.18 


7995 


0.59 


10186 


0.57 


10857 


Vr 


VULCAN/2D 


-0.18 


152.9 


-0.15 


157.6 


-0.1 


161.2 




ZEUS/2D 


0.0 


143.5 


0.26 


156.6 


0.18 


160.9 




VULCAN/2D 


-161.2 


161.2 


-164.0 


164.0 


-164.5 


164.5 




ZEUS/2D 


-146.3 


146.3 


-158.8 


158.8 


-164.3 


164.3 


Br 


VULCAN/2D 


-160.7 


160.7 


-205.8 


205.8 


-214.0 


214.0 




ZETTS/2D 


-127.7 


127.7 


-189.3 


189.. 3 


-208.0 


208.0 




VULCAN/2D 


-0.8 


381.3 


-3.18 


405.9 


-4.5 


417.0 




ZEUS/2D 


-7.8 


331.1 


-6.3 


402.2 


-8.6 


420.0 



TABLE 4 

Energy components in the high-resolution non-rotating implosion test. 



Time 


Code 


Ekin 


Eint 


Emag 


Erot 


J 


Etot 


0. 


VULCAN 


0.0 


4.705 X 10"' 


2.133 X 10^ 


0. 


0. 


4.918 X 10"' 




ZEUS 


0.0 


4.705 X 10"' 


2.117 X 10^ 


0. 


0. 


4.916 X 10"' 


0.03 


VULCAN 


6.575 X 10^ 


4.625 X 10"' 


2.223 X 10^ 


0. 


0. 


4.913 X 10"' 




ZEUS 


6.954 X 10^ 


4.613 X 10" 


2.216 X 10^ 


0. 


0. 


4.905 X 10"* 


0.05 


VULCAN 


6.598 X 10^ 


4.618 X 10"' 


2.270 X 10^ 


0. 


0. 


4.911 X 10"' 




ZEUS 


6.961 X 10^ 


4.605 X 10"' 


2.267 X 10^ 


0. 


0. 


4.902 X 10"' 


0.07 


VULCAN 


1.598 X 10^ 


4.658 X 10"' 


2.353 X 10^ 


0. 


0. 


4.910 X 10"' 




ZEUS 


1.881 X 10^ 


4.645 X 10" 


2.364 X 103 


0. 


0. 


4.900 X 10" 



0.10 VULCAN 6.247 x 10^ 4.624 x 10"' 2.215 x 10-"' 0. 0. 4.908 x 10"' 

ZEUS 6.469 X 10^ 4.611 x 10"' 2.216 x 10^ 0. 0. 4.898 x 10"' 
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TABLE 5 

Energy components in the high-resolution rotating implosion test. 



Time 


Code 


Ekin 


Eint 


Emag 


Erot 


J 


Etot 


u. 


\TJJT P A 1\T 




/I 70 V 1 n4 
1. / \ >o X L\ > 


Z. i-).) X lU 


i .u / n X L\} 








ZEUS 


0.0 


1.705 X 10^ 


2.117 X 10-'' 


1.07(> X 10- 


«7.()8 


1.927 X 10* 


0.03 


VULCAN 


6.199 X 102 


4.626 X 10^ 


2.222 X 10^ 


1.357 X 10^ 


87.08 


4.923 X 10* 




ZEUS 


6.546 X 10^ 


4.615 X 10^ 


2.214 X 10^ 


1.362 X 10^ 


87.70 


4.916 X 10* 


0.05 


VULCAN 


5.861 X 10^ 


4.617 X 10* 


2.255 X 10^ 


1.988 X 10^ 


87.08 


4.921 X 10* 




ZEUS 


6.172 X 10^ 


4.606 X 10* 


2.253 X 10^ 


2.021 X 10^ 


88.68 


4.913 X 10* 


0.07 


VULCAN 


1.693 X 10^ 


4.645 X 10* 


2.323 X 10^ 


2.610 X 10^ 


87.08 


4.920 X 10* 




ZEUS 


1.869 X 10^ 


4.632 X 10* 


2.332 X 10^ 


2.733 X 10^ 


90.16 


4.912 X 10* 



0.10 VULCAN 5.811 x 10^ 4.623 x 10* 2.217 x 10^ 1.600 x 10^ 87.08 4.919 x 10* 
ZEUS 6.037 X 10^ 4.611 x 10* 2.213 x 10^ 1.737 x 10^ 93.10 4.910 x 10* 
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Fic;. 2. — Template of the advection step: the velocity vector defines two donor-type diagonals (a,d) and two acceptor-type diagonals 
(b,c), according to the upwind stability criterion. 
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Fig. 3. — Illustration of the advected Bj,— field for the pure advection test of a transverse magnetic pulse at t = 50, employing 500 zones, 
using the ZEUS code (green: Donor cell scheme; red: van Leer scheme) and VULCAN/2D (solid black). For comparison, we overplot the 
"exact" advection (translation) of the initial pulse (broken black). 
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Fig. 4. — Illustration of the density (top), x-velocity (middle), and S^-field (bottom) for the Brio-Wu test at i = 10, employing 512 
zones, using VULCAN/2D (black) and ZEUS/2D (red). 




Fig. 5. — Left: Illustration of the Li error norm computed on the density for the Brio-Wu test with VULCAN (black) and ZEUS/2D 
(red) as a function of the logarithm of the number of zones employed, i.e. N =64, 128, 256, 512, 1024, and 2048. 
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-1.0 -0.5 0.0 D.5 1.0 -1.0 -0.5 0.0 D.5 1.0 

X X 



-213 -128 -42 42 128 313 83 167 ' £51 335 419 




Fig. 7. — Colormaps of the high-resolution explosion test at t=0.035: density (top left), pressure (top right), Br (bottom left) and 
i?z (bottom right). Each panel shows the ZEUS/2D (left) and VULCAN/2D (right) results. Velocity vectors are overplotted in white 
in the top-left box, with a maximum length set to 10% of the width of the panel, and a maximum magnitude of 165 (all quantities are 
dimensionless). Magnetic field lines are overplotted in black in the top-right box, with footpoints at z = — 0.95rniax equally spaced along 
X every O.lrmax. 
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Fig. 8. — Colormaps of the high-resolution non-rotating implosion test at t=0.05, just before the shock bounces at the inner boundary: 
density (top left), pressure (top right), Br (bottom left) and (bottom right). Each panel shows the ZEUS/2D (left) and VULCAN/2D 
(right) results. Velocity vectors are overplotted in white in the top-left box, with a maximum length set to 10% of the width of the panel, 
and a maximum magnitude of 15 (all quantities are dimensionless). Magnetic field lines are overplotted in black in the top-right box. 
Note the top-bottom symmetry of the z-component of the magnetic field, and the anti-symmetry of the r-component, best rendered by the 
morphology of the poloidal field lines (top right panel). 
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Same as Fig.fS] but at time t = 0.07, corresponding to a time after bounce of ~0.01 seconds. 




-50.0 -20.0 10.0 40.0 70.0 100.0 -150 -90 -30 * 30 90 150 




Fig. 10. — Colormaps of the rotating implosion test at t=0.05 and t=0.07: The angular velocity Q (left column), the toroidal magnetic 
field (right column), at times t = 0.05 (before bounce; top row) and t = 0.07 (after bounce; bottom row). Each panel shows the 
ZEUS/2D (left) and VULCAN/2D (right) results. Again, note the top-bottom symmetry of the z-component of the magnetic field and the 
anti-symmetry of the r-component, best rendered by the morphology of the poloidal field lines (right panels). 
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Fig. 11. — Colormaps of the high-resolution non-rotating implosion tests with dominant initial magnetic fields {Bz = 120, t = 0.06, top 
row; Bz = 300, t = 0.04, bottom row). Shown are the density (left column) and the quantity /3piasma = ^'/C^^/S'"") (right column). Each 
panel shows results for the ZEUS/2D (left) and VULCAN/2D (right). Again, note the top-bottom symmetry of the z-component of the 
magnetic field and the anti-symmetry of the r-component, best rendered by the morphology of the poloidal field lines (right panels). 



